Estimation of classical parameters via continuous probing of complementary quantum 

observables 
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We discuss how continuous probing of a quantum system allows estimation of unknown classical 
parameters embodied in the Hamiltonian of the system. We generalize the stochastic master equa- 
tion associated with continuous observation processes to a Bayesian filter equation for the probability 
distribution of the desired parameters, and we illustrate its application by estimating the direction 
of a magnetic field. In our example, the field causes a ground state spin precession in a two-level 
atom which is detected by the polarization rotation of off-resonant optical probes, interacting with 
the atomic spin components. 
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I. INTRODUCTION 

High precision metrology with quantum systems is a 
research field of high current activity. Through the quan- 
tization of their energy levels, elementary quantum sys- 
tems provide fundamental time and frequency standards 
and, due to the highly developed means for preparation, 
control, and detection of these systems, they serve as ex- 
cellent probes of various perturbations such as applied 
electric and magnetic fields, and incrtial effects associ- 
ated with rotation, acceleration, and classical or rela- 
tivistic gravitational effects. 

The theoretical research proceeds along different di- 
rections according to the different measurement schemes. 
Thus, for experiments where a quantum system is subject 
to a perturbation for a given short duration of time, the 
search for the initial quantum states on which different 
values of the perturbation leads to the most distinguish- 
able outcomes has promoted the use of concepts such as 
the Fisher information and the Cramer- Rao bound [![, 
and has identified squeezed states, Schrodinger cat-like 
states and generalizations hereof as useful resource states 
in metrology. Along a different path, measurements that 
occur continuously in time, such as continuous wave laser 
spectroscopy, are made subject to analyses, that serve to 
exhaust the information about the desired parameters 
from the entire sequence of measurement data. While a 
simple relationship between the average signal, e.g., an 
absorption profile, and an unknown physical parameter 
provides a relatively straightforward method for estima- 
tion, the systematic extraction of a reliable error-bar on 
the estimate is more challenging [HI]- 

In this work we present such an analysis for the non- 
trivial case of the continuous probing of a single quan- 



tum system. We detail a Bayesian analysis, which treats 
our description of unknown parameter values and the un- 
known state of the quantum system on an equal footing 
and where any data received serves to update our prior 
estimate of the parameters of interest. This idea has been 
previously applied to the case of quantum-non-dcmolition 
(QND) probing of quantum systems [3, H|, where it is 
largely equivalent to a Kalman filter [3(. The formalism 
presented here, however, is more general, and new physi- 
cal features appear, when not only a single QND variable 
is being detected. 



The paper is organized as follows: In Sec. II we present 
our general quantum filter equation and we give an ex- 
plicit derivation for our specific physical model system. 
In Sec. Ill we introduce the unknown classical parame- 
ters and we show how their representation as effective 
incoherent quantum degrees of freedom augments the 
quantum filtering equation to automatically provide a 
Bayesian update formula for the unknown classical pa- 
rameters. In Sec. IV we present numerical simulations 
and we show how the fluctuating measurement signals 
of optical probes interacting with an atomic spin gradu- 
ally filters the probability distribution of a magnetic field 
with known strength, but unknown direction, and how it 
ultimately reveals this direction with good precision. In 
Sec. V we conclude and discuss the outlook and perspec- 
tives of our work. Details about the derivation of the 
augmented stochastic master equation (SME) and on its 
numerical solution are provided in the appendix. 
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II. OPTICAL PROBING OF A SINGLE ATOMIC 

SPIN 

A. General quantum filtering equation 

A quantum system with Hamiltonian H subject to 
Markovian damping, described by Lindblad operators Oj 
and rates Tj, is described by a reduced system density 
matrix g s , which obeys the master equation 



dg s = \^[g s ,H] +Y i r j T>[d j ]g a dt, (1) 



where the operator T> is defined as 



-Dim = HiP 



PfQl + QlPf 



(2) 



Interaction with continuous quantized probe fields 
cause entanglement of the system with the fields which, 
if the field degrees of freedom are subsequently traced 
out, is described by inclusion of further master equation 
terms T>[M. n ]g with Lindblad operators M. n and effective 
interaction parameters M n . The entanglement of the sys- 
tem and the fields, however, implies that measurements 
of the probe field variables after the interaction lead to a 
back-action on the state of the quantum system. 

Subject to the quantum back-action due to continuous 
amplitude measurements on the probe field, the reduced 
density matrix of the quantum system obeys the follow- 
ing stochastic master equation 



dg s = \ ^[q s ,H) + I'/PC', "' + X] M n V[M n )g\ dt 
+J2 V*lnM n H [M n ] g s dW n (t) , (3) 

n 

where the operator T-L is defined as @ 



H[fW = fg s + g a p-{f + P)@ 



(4) 



and where dW n are infinitesimal Wiener processes ac- 
counting for the noisy contribution to the field amplitude 
measurement 



dY°{t)= Vn 



Ml)dt + dW n {t). (5) 



In the equations ((3]) and (j&J) the parameters r\ n account 
for the detector efficiencies and transmission losses of the 
probe beams between the system and the detector. If 
Vn = 0, the corresponding probe field merely contributes 
extra damping and dccohercncc to the system, while if 
rj n = 1, and in the absence of other damping or inef- 
fective probing terms, the system may be described by a 



stochastic wave function rather than a density matrix 0- 
[Toj . The random Wiener noise increments dW n consti- 
tute the so-called innovation processes, i.e., they are the 
difference between the experimentally observed signals 
dY®(t) and their quantum-mechanical expectation val- 
ues r\ n \/ M n (M. n + Mj^dt determined from the current 
quantum state of the system. This difference amounts to 
the shot-noise in field amplitude measurements by homo- 
dyne detection, and the formalism can also be adapted 
to treat, e.g., photon counting measurements, where the 
innovation process is described by an (infinitesimal) Pois- 
son process. In either case, the measurement back-action 
has only infinitesimal influence in every small time step 
while for certain measurements accumulated over time 
it ultimately causes the collapse on a random eigenstate 
normally attributed to the von Neuman projection pos- 
tulate. 

We have presented the quantum filtering equation in 
a general form following a Markovian and perturbative 
treatment of the interaction of the system with its envi- 
ronment. For every concrete physical example one has to 
validate such a treatment and to explicitly analyze the 
appropriate interactions and field measurement schemes 
to obtain the operators and parameters in Eq. ([3]). In 
the next subsection we will recall such a derivation for 
the explicit case of an atom with a degenerate ground 
state that interacts with off-resonant probe laser fields. 



B. Dynamics of a two-state atom and an 
off-resonant laser field 



The detection of optical phase shifts and polarization 
rotation is at the heart of spin squeezing and quantum en- 
tanglement schemes involving the collective spin degrees 
of freedom associated with large atomic ensembles [Til — 
[T^ |. and it also constitutes the basis for atomic mag- 
netometers which today explore the quantum limits of 
resolution (l5l - [l7j . The spins in large atomic ensembles 
are well approximated as harmonic oscillator degrees of 
freedom and Gaussian approximations and classical fil- 
tering theory apply [H, O, [l|| , while single atoms must 
be described by the full density matrix formalism, that 
we will address in the following. 

We consider an atomic quantum system with degen- 
erate ground states \g m ) and excited states |e m ), where 
m = ±1/2 denotes the azimuthal quantum number with 
respect to the quantization axis z. The atom interacts 
through its magnetic moment /2 with a classical magnetic 
field B = (B x , B yi B z ), which drives a Larmor precession 
of the ground state atomic spin. The atom is coupled to 
an off-resonant linearly polarized laser field. This field, 
which is linearly polarized along the x-axis, can be de- 
composed into two circular components with annihilation 
operators 6± = (T^x + ia y )/\2. Due to the dipole selec- 
tion rules, the two circularly polarized field components 
couple individually to the two different ground state pop- 
ulations. 
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During off- resonant probing (i.e., A ^> g), the atomic 
excited state can be adiabatically eliminated, and the 
quantized field-atom Hamiltonian is given by [20 ]: 



a\at\g(/ 2 ){gt/2\ +p-B. 



(6) 



e=±i 



In Eq. ([6]), A = u>l — is the laser atom detuning, 
g = d-Eo/h, where d is the atomic electric dipole moment 
and |J3o| = \J hhj l / {V eo) is the electric field per photon 
in the quantization volume V, and eo is the (electric) 
vacuum permeability. 

As as result of Eq. ((6J, the two circularly polarized 
field components experience phase shifts that depend on 
the atomic occupation of the two ground states. The 
resulting phase difference between the two field compo- 
nents implies a (Faraday) rotation of the field polariza- 
tion, which is proportional to the population difference 
between the ground states (#±1/2)- 

Because of the strong linearly polarized probe field 
with photon number iV p h S> 1, the Stokes operators of 
the field can be written as 



J, 



a x a x 
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(&„ +a+) = y/N^y, 



{fly - a\) = V^VphPy, (7) 



defining the canonical conjugate operators y and p y . The 
polarization rotation of the field is conveniently measured 
by subtracting the intensities of polarization components 
linearly polarized at ±45 degrees with respect to the in- 
cident field, and when this quantity is expressed in terms 
of the Stokes obscrvables we recover an expression pro- 
portional to the operator y. 

By writing /2 = /act, where a = (a x ,ay,a z ) are the 
Pauli matrices, and using the definitions in Eq. ([7]), the 
total Hamiltonian can be written 



segment, in turn, is described as a single harmonic oscil- 
lator mode described by the operators in Eq. ([7]). In the 
time interval t the atom interacts with iV p h = <frr pho- 
tons, where <E> is the photon flux, and we assume that the 
interaction is sufficiently weak, that the dynamics can be 
well described by the coarse-grained propagator 



U = e 



where 



hg 2 r 



phi 



/IT, 



(9) 



(10) 



and Hb = { n Bi n Bi n B) 18 the unit vector pointing in 
the magnetic field direction. We have introduced &a B — 
a ■ Hb and, in the part describing the atom-light inter- 
action, we have introduced = a ■ m, for the Pauli 
operator along an arbitrary unit vector direction fn. To 
this aim we assume that we can probe the atomic ground 
state spin along any such direction with probe beams of 
appropriate polarization and propagation direction, or 
by application of a unitary spin rotation prior to prob- 
ing of the z-component with a fixed beam set-up. Since 
iV p h oc t and g oc t" 1 / 2 (through the volume V = Arc), 
the dimensionless coupling constant n T is proportional to 
r 1 / 2 , while /i r is linear in r, and in Eq. we have thus 
neglected terms of order r 3 / 2 and higher. 

The incident laser beam is in a coherent state of lin- 
early polarized light described by a Gaussian wave func- 
tion 7r~ 1 / 4 e _p «^ 2 in the argument p y , associated with 
the observable p y introduced in Eq. ((7J). The joint 
state of an atomic ground superposition state \ip 8 (t)) = 
J2i=±i c i/2\9i/i) an d the incident y-polarization compo- 
nent of the quantized probe field can thus be written in 
the product basis \p y ,gi/ 2 ) = \p y ) ® \gtj2), 

I*p.(*)> = 4j4 E c m I dPye- pl/2 \Py,9i/2)- (11) 

Under the action of (j9]), this state evolves into the entan- 
gled state 



(8) 



a—x,y,z 



A term ^(a x a x + d y d y ) has been omitted from Eq. as 
it gives rise to a common phase shift, but no polarization 
rotation of the probe laser field. 

To properly account for the entanglement created be- 
tween the atom and the continuous probe field, we treat 
the laser beam as a sequence of segments of length 
L = ct, area A = V/L and volume V, each initially 
prepared in a coherent state before the interaction. The 
continuous interaction between the atom and the light 
beam can then be represented as a sequence of interac- 
tions of the atom with one segment after the other. Each 



|#pb(* + t)} = e 



= IT/I Yl I d Pv c 'i/2{Py) e p2y/2 \Py, 



9t/2)- 

(12) 



where, to first order in t (second order in y/r), the new 
expansion coefficients c' e , 2 (p y ) are: 



"1/2 = C <72 



h 



flr\B\ 



t-C-112 



-j-Pyfax - iimy) + ^Y^-{n B - i£n v B ) 



(13) 
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C. A quantum filtering equation for the two-state 
atom 

To describe the back-action due to the field measure- 
ment, it is convenient to transform the entangled state 
(|T2|) to the y rather than p y representation of the field, 
using the relation 



\Py) 



/2tt 



dye 



-iypy 



\y) 



(14) 



By replacing y with the expression (|19[) for y D in 
Eq. (JTTJ) , and by expanding the expressions for the state 
amplitudes to lowest order in t, we obtain, in the contin- 
uous limit, the quantum filtering equation for the state 
of the atomic system 



dg s 



i— [<rg, g s ]dt + MT>[arh]g s dt 



liji 
h 

MH[a A }Q s dW(t). 



(20) 



and the state (fl2l) can be rewritten as: 



Here the interaction parameter M is given explicitly by 



3 y /2 \y,gi/2) 



(15) 



with the new coefficient Q/ 2 (y) given by 



Q/2 = lC-£/2 

1 /K^ 2 



ce/2 



i£\n 



, Mr|B| 



vm z —y 



(16) 



A measurement of the light probe observable y with 
outcome y D projects the state of the system (fl~5|) onto 
the state component with that definite value, i.e., the 
atomic part of the system becomes 



(17) 



e=±i 



where the explicit dependence of q/ 2 on the measure- 
ment outcome y D causes an (infinitesimal) transfer of 
amplitude among the two atomic states. 

Given the quantum state (jT5j) , the probability to mea- 
sure a given value y D is 



V(y D ) = £ \c e/2 (y 

v n e=±i 



12 ^ 7T -l/2 e ^(y D -y ) 2 



(18) 



where yo = ^-(o'm)- This explicitly shows how the op- 
tical probing yields a signal proportional to the desired 
mean value (<r,n), and it allows us to to model the mea- 
surement outcome y D as a stochastic variable: 



d k t AW 

y =y( (J ™/" 



2r 



(19) 



where AW is a (finite) Gaussian Wiener increment with 
mean zero and variance r. 



M = 



9 4 t 2 



4(w L - uj A y 



(21) 



and an infinitesimal Wiener process models the noise in 
the detected signal, dY D (t) = 2y/M(a m )dt + dW. Thus, 
the system evolves according to a stochastic equation of 
the same form as the standard quantum filter equation 
©• 

We note that several probe fields may be used to probe 
different spin components, and, to lowest order in r, their 
effects on the quantum state commute. They may hence 
be included as separate terms with independent Wiener 
processes dW n . 

The modifications in Eq. due to finite detector effi- 
ciency can also be understood from first principles in the 
model system, and lead to similar correction factors in 

Eq. m. 



III. CONDITIONAL DYNAMICS AND 
ESTIMATION OF A CLASSICAL PARAMETER 



We are interested in the use of quantum systems to 
estimate a classical physical parameter or a set of param- 
eters 7. Here, in order to keep the discussion as general 
as possible, we treat the unknown parameter 7 as a vec- 
tor quantity to indicate that it may be a set of parame- 
ters such as a damping rates, energy shifts, and coupling 
strengths, or, as in our example below, the directional 
components of a vector magnetic field. The experiment 
is sensitive to the value of these parameters, e.g., if they 
are coefficients in the Hamiltonian H = ^(7) acting on 
the system, and if this dependence results in a change of 
the observables probed in the experiment. 

We describe the quantum dynamics of the combined 
system with 7 belonging to a finite set of values V 7 = 
tf k :k = l,...,N}. 

For an observer who knows the true value of 7 = 7fe > 
the system is described by our original reduced system 
stochastic master equation ^ with H = ff(7fc ), 
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dg s = Ug s , H(n )}dt + £ T 3 V[6 }el&t 



A I 



+ ^ M n V[M n ]g s dt + ^r] n M n %[MnWodW n {t), 

n=l 

(22) 

conditioned by the measurement signals 



dY°(t) = ^dW n (t) 



(23) 



where (M n + M ] n ) Q = Tr s {(7W„ + MUqq}, and where 
dW n (t) arc standard Wiener processes. Here Trs(-) is 
the trace on the system Hilbert space. 

In the following we will assume that we probe all 
three spin components {cr x , a y , a z ) with measurement 
strengths (Mi,M 2 ,M 3 ) and efficiencies (771 , 772 , 773 ) . By 
introducing the three-dimensional real Bloch vector, r = 
Tr(gQ<j), the stochastic master equation (f2"2"j) can be 
rewritten in the Bloch vector representation, 



dr = 2 I (b k „ x r ) - (ai + a 2 + a 3 )f + 
+ (l-r 2 dfi)-fx (rxdfi)}, 



(24) 



where we have passed to dimensionless units by 
performing the replacement t — > Alt with M = 
max{Mi, M2, M3}. Besides this, we have defined the vec- 
tors a, with components a n = M n /M, ff = (771,772,773), 
bk = H-BBk/ (hM), and dfi, with components dtt n = 
y/ri n a n dW n . The symbol * in Eq. ([24]) indicates the 
pointwise product, (p * q) n = p n q n - 

Equatio n (1241) was derived and analyzed in detail in 
Refs. [U |22|. for the special case of identical probing 
strengths and efficiencies for all three orthogonal spin 
directions. As shown in that work, in the absence of 
a magnetic field, the Bloch vector is driven towards a 
steady state radial (purity) distribution and an isotropic 
angular distribution within the Bloch sphere. The mag- 
netic field, in turn, provides a torque for the atomic spin 
vector, and the resulting spin precession, in the plane 
perpendicular to the field, reveals itself in a modulation 
of the polarization rotation measurements. 



A. Augmented quantum filter equation 

In our formulation of the estimation problem we will 
treat the classical parameter 7 as unknown with an as- 
signed probability distribution P{^). The measurements 
then cause an update of the probability distribution, 
which is governed by Bayes rule for conditional prob- 
abilities. Until the actual value of 7 is known, we thus 
have to treat each candidate value with a probability fac- 
tor, and for each possible value of 7, the corresponding 



state of the quantum system evolves under the quantum 
filtering equation with the corresponding dependence of 

7- 

To this end it is convenient 0, IH [23Tj25j to consider 
the augmented Hilbert space fj = fj s <£) -f) 7 , where fj s 
and f) 7 refer to the quantum system Hilbert space and 
the space of classical states for the variable % respec- 
tively. The latter space describes states with definite val- 
ues of the parameters, and superposition states are not 
populated. The quantum mechanical notation, however, 
still applies and, e.g., describes a probability distribu- 
tion for a set of values 7*fc as a diagonal density matrix 
g> 7 = Y^k Pk\lk)(lk\- The classical variables 7 are equiv- 
alent to quantum non-demolition (QND) variables of an 
ancillary quantum system that interacts with our probe 
system, and for which the Baycsian probability update 
is fully equivalent to the quantum measurement back- 
action. When we incorporate the parameters 7 in this 
way we can directly apply the filtering equation on the 
augmented space. 

The observer who does not know the value of 7 de- 
scribes the combined quantum and classical system by 
the augmented density matrix 



JV 



(25) 



k=l 



where g B k is the normalized system state density matrix 
associated with the specific value 7 = 7fc- 

The combined system evolves according to the quan- 
tum filter equation 



dg= -[g,H}+J2^MOj}g + ^2M n V[M n }g\ dt 



n=l 



+Y J ^M~nU[M n ]g (dY°(t) - r} n ^fM^{M n + Ml) E dt) , 

(26) 

where the operator H is defined as: 



n[f}g = fg + gP-(f + P)Eg. 



(27) 



Here we use the notation (/) e = Tr(/g) = 

J2k=i Pk^s{fg%) to explicitly recall that the expecta- 
tion value of the signal should be determined by the full 
augmented quantum state, equivalent to a weighted av- 
erage over the ensemble of states of the quantum system 
governed by the different values of 7. 

If the dependence on 7 only enters via the Hamiltonian 
^(7), the different terms in (|26|) are implemented as the 
following product operators on = S) s <8> ,f) 7 



/v 



(28) 



fc=i 



6 



N 



fe=i 

N 

M n = \lk)(lk\ ® M n - 



(29) 



where each value of 7& is associated with an unnormalizcd 
p s k . Equation (J32J) is of course equivalent to Eq. (|25|) 
with the normalized system density matrix gf. and the 
probability distribution P k . 



k=l 



B. Detection signal properties 

The stochastic process appearing in Eq.(|26|). 

7^dV n {t) := dYf(t) - Vnx/K(M n + Ml) E dt. (30) 



is not a standard Wiener process. This is because we sub- 
tract from the measured signal a weighted average based 
on our probabilistic description of 7, while the measured 
photocurrent dY® in the experiment is governed by the 
actual value of the unknown parameter 7 = . 

Equation (|23|) characterizes the properties of such a 
realistic detection record, and, when inserted in Eq. (f2l>|). 
we find that the stochastic process in Eq. (|30|) can be 
rewritten as 



As shown in detail in the appendix IV A[ Eq. (f2l))) leads 
to two separate equations for gt and 



dgl = \ [Q%,H(%)]dt + J2 rMOMM 

3 

M 

+ Y, M n V[M n }g\dt + ^M~ n n[M n ]g s k dV^(t), 

n=l 

(33) 



with dv]i k \t) = dY°(t) - rin^M^(Mn + Ml) k dt, and 



dP k = Tr s (dpi) 



A I 



dV n =dW n -^r] n M n [(M n +M\) E - (M„ + M n ) \dt. = P^Vm^ [{M n + Ml) k - {M n +Ml) E 

(31) n=l 



dV n (t). 
(34) 



Hence, dV n (t) is a stochastic Gaussian process with vari- 
ance dt, and with (statistical) mean value ((dV n (t))) = 
y/Vn M n ((M n + M\ L ) -{M n + M\ l ) E )dt, reflecting pre- 
cisely the difference in the expectation value assumed by 
the weighted average over different j k and by the correct 
value 7fc . 

Equation (|26|) permits a full simulation of the detection 
process and provides a time dependent solution of the 
form 



JY 



= 5Zl7fe)(7fc|®/3|, 



(32) 



k=l 



In both equations the photocurrent dY^(i), observed or 
simulated according to Eq. (|22|) . appears, and Eq. ([34|) 
agrees with the expression given in Ref. |2J]. We note, 
however, that the stochastic term in our Eq. (f3"3")l con- 
tains the expectation value (M n + Mj)k corresponding 
to the parameter value j k , while in Ref. 24| the ensemble 
average (M n + M n )E has been used. 

By inserting the expression (|3T|) for dV n in Eq. (|34|) . 
we see that the change of dP k due to the measurements 
is given by 



(M n + Ml)k - (M n + MDe) dW n (t) + v^A4 (M n + MDo - (M n + M n ) E ) dt 



(35) 



r 



This equation has a natural interpretation: For pa- 
rameter values jk where the expected mean current 
oc (M n +Ml t )k differs in the same (opposite) direction 
from the ensemble mean as the one expected for the ac- 
tual value (M n + M n )o, the two parentheses will typi- 
cally have the same sign, and P k will increase (decrease). 
Due to the random contribution dW n (t), however, the 
probabilities will show fluctuations, and their increase 



(decrease) with time will appear only as an average trend, 
leading, in particular, to a typically increasing value for 
the probability of the correct value P ko . 
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IV. VECTOR MAGNETOMETRY WITH A 
TWO-LEVEL ATOM 

We now apply the formalism to the two-level atom cou- 
pled to a magnetic field B with known magnitude \B\, 
but unknown direction in space. Since such a field will 
cause a spin precession around the magnetic field axis, 
we expect that optical probing of a single spin compo- 
nent will not be sensitive to the magnetic field projection 
along the spin direction probed, while the simultaneous 
probing of all three spin components is bound to reveal 
any motion of the mean spin vector due to the magnetic 
precession. 

The acquisition of data from a real or a simulated ex- 
periment cause a continuous update of the probability 
distribution Pk for the magnetic field, represented in the 
following by the dimensionless vector, bk = (J>b Bk / (flM) . 
The angular measure, cos0(i) = \b u \~ 2 J^k Pk(t) (bk ■ bu), 
quantifies the scattering of the magnetic field directions 
bk inferred from a single experimental run around the 
actual value b u . By carrying out a large number of simu- 
lations, we thus quantify the average performance of the 
method by the (average) scalar product 

«cos0»(i) = \b u \- 2 ({ p k (*) (h ■ b u ))) (36) 
k 

as a function of the measurement time t. 



A. Direction of a magnetic field along a given axis 

Following Ref. |23 |. we have first investigated the case 
in which the initial state of the atom is represented by 
the Bloch vector f = (0, 1, 0) with M h2 = but M 3 ^ 
(i.e., we probe only the component of the spin along z and 
M = M3), and where the magnetic field has a known 
strength while it has equal prior probabilities to point 
in the positive and in the negative a;-axis directions. In 
Fig. [T] we show the behavior of ((cos #))(£) as function of 
time in the two cases of a weak |6„| = 0.1 (lower, black 
curve) and a strong \b u \ = 1.5 field (upper, red curve). 
The curves are obtained by averaging over 104 simulated 
detection records. The stronger the field amplitude, the 
better is the directional estimate, but, as also observed 
in Ref. [24| . the quantum filter does not unambiguously 
identify the direction b u of the unknown field. 

This situation changes when all the three spin compo- 
nents are detected. As illustrated in Fig. [3J the accumu- 
lation of results from all three detectors lead, especially 
for high strength of the field, to a final probability distri- 
bution which is well converged to the correct b u . In these 
numerical experiments 104 trials have been performed in 
order to accumulate sufficient data for the statistical av- 
erages within a reasonable computational time. A small 
fraction (~ 1%) of the simulated trajectories are unsta- 
ble in the case of the three detectors and they have been 
rejected from the statistical averages. 
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Figure 1: (Color online). Time evolution of ((cos #}}(£) for 
b u = (0.1,0,0) (lower, black line) and b u = (1.5, 0, 0) (upper, 
red line) for the qubit example discussed in Ref. [24 ]. 

There is a number of competing effects that may ex- 
plain the dependence on the quality of our estimate on 
the field strength and the number of spin components 
probed: Measuring the z-component of a spin may lead 
to a (Zeno-effect) suppression of, and, hence, insensi- 
tivity to, slow precession of the spin around the ir-axis, 
while even for strong fields, the measurement of a single 
spin component does not allow discrimination between 
left and right circular precession around the x-axis. The 
probing of several components on the one hand allow 
discrimination between left and right circular precession, 
and, on the other hand prevents (Zeno-) locking of the 
system to eigenstates of the probed quantities as they do 
not commute and do not have common eigenstates. 



B. Estimation of a spherically random direction of 
a magnetic field 

Finally, we have analyzed the scenario in which the 
magnetic field has an isotropic prior probability distribu- 
tion, represented by an ensemble Vb with TV = 98 direc- 
tions on the unit sphere, as illustrated in Fig. [3] Here, 
we assume the stronger field \b u \ = 1.5, and we assume 
equal strength probing of all three spin components. 

The initial condition for the probabilities with all 
Pk = l/N is illustrated by the (green) sphere displayed 
in Fig. [3] at time Mt = 0. We study the convergence of 
the probability distribution as a function of the probing 
time, and for long times (MT = 15), the filter converges 
well to a single direction. We note that the spheres are 
cut on the top because in the simulation we have consid- 
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Figure 2: (Color online). Same setup as in Fig. [T] but with 
three detectors (M = Mi = Ah = M 3 ). 



ered an interval 9 = (0, tt) equally spaced with Ng = 7 
grid points and an interval 4> = [0, 2tt) with Afy = 14 grid 
points for the azimuthal angle. 

The equal strength probing of the three cartesian spin 
components cr x , a y , a z is equivalent [2ll l24j to prob- 
ing of any other cartesian set, including for example one 
parallel component and two perpendicular components 
to the applied magnetic field, and we find that the moni- 
toring of all three spin components lead to unambiguous 
identification of the direction of the applied field. 

The spin vector may, inadvertently, align, parallel or 
anti-parallel with the applied magnetic field axis, but the 
random back-action of the probing of the non-commuting 
orthogonal spin observables will kick the system away 
from these directions and give rise to renewed observable 
precession. 



ics. In this mapping, we model the classical parameters 
as QND observables of auxiliary quantum systems, and 
their classical probability distribution thus coincides with 
the conventional reduced density matrix elements for a 
quantum system. Since the quantum state description 
cannot be completed by further knowledge in the form of 
hidden variables, our formulation of the parameter esti- 
mation problem, indeed, provides the tightest and most 
precise probability distribution for the variable of inter- 
est conditioned on the measurement outcome and on the 
prior probability distribution. 

The discretisation of the parameter space and solution 
of a quantum system master equation associated with 
each potential parameter value naturally puts limit on 
the precision of the method and the number of variables 
that can be realistically determined. A natural next step 
would be to apply methods that gradually refine the pa- 
rameter space around the most likely values and sup- 
press the most unlikely ones from the calculation. Such 
weighted stochastic differential equations are known in 
statistics (26j . and we imagine that they may be used to 
decide objective means to suppress and to breed new pa- 
rameter values without enlarging the memory and com- 
putational demands of the method. Another interesting 
approach, put forward in Ref. [27j . involves projection of 
the complete system on a non-linear lower-dimensional 
manifold on which the integration of the stochastic differ- 
ential equations of motion is faster. Alternatively, max- 
imum likelihood methods and random searches through 
the parameter space, e.g., by Markov Chain Monte Carlo 
methods [28j . may be effectively applied to even very 
large search spaces. We imagine that our simulations 
may serve as useful reference data for testing such alter- 
native estimation techniques. 



V. CONCLUSION 

In this work, we have demonstrated a Baycsian filter 
for classical parameters, which affect the dynamics of a 
quantum system. Previous studies along the same lines 
have focused on quantum non-demolition measurements 
of typically a single variable, but as shown by our anal- 
ysis, a non-QND setting may be analyzed by the very 
same assumptions and methods. Non-QND probing may 
have specific advantages and provide more decisive re- 
sults, when the parameters affect different observables of 
the quantum probe, as illustrated explicitly by our nu- 
merical simulations. 

The Baycsian filter is derived from a standard quan- 
tum filter formulation of conditioned quantum dynam- 
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Figure 3: (Color online). Upper panels: Time evolution of the probability distribution Pk for an ensemble Vb of 98 elements 
initially equally distributed over the unit sphere. The unknown magnetic field that has to be determied is the strong field 
b u — (f.5,0, 0), and its direction is displayed by the pink dot on the sphere at time Mt — 15. Lower panel: Evolution of 
{(cosd))(t) for the single shot experiment simulated in the upper panels. 



Appendix 

A. Quantum filtering equations for the parameter 
estimation problem 



In order to derive the equations (|33|) and (|34[) we first 
need the SME for p|. which is given by 
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dpi 
Pk 



A I 



A I 



-[Q%,H(j k )} + E M nV[M n ]gl ] dt + Y^VnXiMnQl + - (M n + M\ ) e qI) dV n (t), 

n=l 

(37) 



where we note that dp s k = d((7fc|(5|7fc)). Given this, the 
equation of motion for the probability P k is easily ob- 
tained by computing the trace of Eq. (I37|) . which provides 
precisely Eq. (f34|) . Now, since g s k = p%/P k we have: 



Bloch vector, and propagate the collection of Bloch vec- 
tors as subject to the noisy detection signals - governed 
by Eq. ([23"]) . Using the same notation as in Eq. 
these equations have the form 



dgl = ■ dp| + pi ■ d (i-) + dpi ■ d (±- ) 
where (classical Ito formula [29| ) 



Thus, we have 



-^dP k + ^{dP k f. 

k k 



(39) 



| (h x f fc ) - (ai + «2 + a3)rfe + a * r fe 



dr fe 



+ (1 - 4 ,fc) (4(1 - rg) - r k x (f fc x 7fc)jj dt 
+ dti(l - r 2 k ) - T? k x (r k x dft)} , (43) 

Instead, the probabilities represented as the column vec- 
tor P = (Pi, P2, . . . , Pn) T obey the following matrix 
equation 



(dP fe ) 5 



pi 

™ n—1 



]T r?„M„ (<M„ + M n ) k - (M n + Ml) E ) dt, 

(40) 



where we used the fact that dW n (t)dt = and 
dW n (t)dW n > (t) = 8 n>n idt. Consequently, we obtain 



A I 



d (p~) = P~^ { ? ?« A/ » + M «)k - (M n + MDe 



n=l 



X ((M n + ^)fc - (M n + M n ) k ^) dt - ^fr) n M n 
x ((M n +Ml) k - (M n + M n ) E ) dW n (t)} , 



(41) 



and therefore 



(1 \ 
7T =E VnM n ((M n + Ml) B - (M n + Ml) 
k ' n=l 

X (M n Ql + Q%Mi - (M n + M n ) E Qh) dt. 

(42) 



Putting all together into Eq. (|38|) and by using the defi- 
nition (l23l) we derive (l33l). 



B. Numerical simulations 

In analogy to Eq. ([2~4")l . we can represent the (normal- 
ized) density matrices associated with each value j k as a 



dP = 4 jp * (w T ■ C) T - P P T • (u T ■ C) T 



P* 



(CP) T C 



P 



(CP) T ■ (CP) 



dt + Gdn. 

(44) 



Here, C and C are 3 x N matrices and G is an N x 3 
matrix containing the Bloch vector solutions of Eq. (|43|) , 



» 



P T ■ X„) where X 



C n j — T~in,QnC n ,ji and Gj iTl — 2Pj(x 



(») 



» 



) T , and 



\ x l i • i 2 j ' ' ' i A Af 

7/ * a * f ko ■ 

For the numerical simulation of both (j43|) and (|44|) we 
employed an Ito-Euler integrator [29| with a time step At 
ranging from 2x 10~ 7 • M -1 to 10~ 5 • A/ -1 depending on 
the size of the set Vb and on the number of switched off 
detectors. Such a choice enabled us to have an efficient 
integrator, even though some of the quantum trajecto- 
ries might have been unstable. To solve the instability 
problem, we have first tried to apply an implicit Milt- 
sticn method [H, & but since both 63 and © are 
nonlinear, one has to solve numerically at each time step 
(e.g., by means of the Nelder-Mead method (28|) implicit 
equations like f k (t + At) = r k (t) + f(f k (t + At)), where / 
is the r.h.s. of Eq. (|4"3")) plus some additional term due to 
the Miltstien routine. While such a strategy might solve 
the problem, we have numerically observed that such 
an approach is significantly more time consuming than 
the Ito-Euler integrator. Thus, we also applied a deriva- 
tive free order 2.0 weak predictor corrector method [30l |. 
which turns out to be quite efficient in the case of a single 
Wiener noise process, but in the case of three detectors, 
whose generalization is not straightforward, we noticed, 
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as for other predictor-corrector methods, that the insta- 
bility could not be fixed. Hence, we employed the simple 
Ito-Euler integrator with (rather) small time steps (up to 
At = 2 x 1CT 7 • M _1 ). We noticed that with such a sim- 



ple strategy the unstable quantum trajectories could have 
been reduced or even suppressed, but at the expenses of 
a very long numerical computation. 
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